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Mathematical  techniques  are  presented  which  allow  for  analytical  solutions  of  the  catalyst  layer  transport 
and  electrochemical  problem  in  PEM  fuel  cells.  These  techniques  transform  the  volumetric  reaction  terms 
to  boundary  flux  terms,  thereby  eliminating  the  need  for  computational  solving  of  the  catalyst  layer  prob¬ 
lem.  The  result  is  a  semi-analytical  fuel  cell  model— a  computational  model  that  entails  analytical  rather 
than  computational  catalyst  layer  solutions.  This  helps  to  alleviate  the  meshing  difficulties  inherent  in  the 
catalyst  layers  caused  by  large  geometric  aspect  ratios,  and  hence  reduce  the  computational  requirements 
for  fuel  cell  models. 

These  analytical  solutions  are  implemented  in  a  3D  PEM  fuel  cell  model,  and  the  results  of  the  semi- 
analytical  model  match  well  with  the  full  computational  model  in  terms  of  the  polarization  performance 
and  species  concentration  distribution.  In  addition,  these  analytical  solutions  were  able  to  reduce  the 
required  computational  memory  by  a  factor  of  approximately  3,  and  the  computational  time  by  a  factor 
of  approximately  4. 

©  2008  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Fuel  cells  have  generated  considerable  interest  in  recent  years, 
being  heralded  as  the  power  generating  devices  of  the  future.  Pro¬ 
ton  exchange  membrane  (PEM)  fuel  cells  are  prime  candidates  to 
eventually  replace  the  internal  combustion  engine  (ICE)  in  vehicles. 
They  use  hydrogen  as  fuel,  which  reacts  electrochemically  (without 
direct  combustion)  with  oxygen  to  produce  electrical  power.  The 
only  emissions  from  a  fuel  cell  are  water  and  heat,  thereby  making 
them  more  environmentally  friendly  than  ICEs. 

Fuel  cell  modeling  and  simulation  are  critical  aspects  of  the 
development  of  fuel  cell  technology.  This  is  evident  from  the 
numerous  publications  pertaining  to  fuel  cell  modeling  over  the 
past  15  years.  Fuel  cell  modeling  has  evolved  from  ID  to  3D  [1-3], 
from  steady  state  to  transient,  from  single  phase  to  two  phases 
[4],  and  from  straight  channels  to  more  complex  serpentine  and 
inter-digitated  flow  fields  [5-8].  Our  group,  as  well  as  others,  were 
responsible  for  the  development  of  intermediate  to  high  tem¬ 
perature  PEM  fuel  cell  modeling  using  alternative  membranes  to 
Nation®,  e.g.  polybenzimidazole  (PBI)  [9-16]. 


Abbreviations:  CL,  catalyst  layer;  GDE,  governing  differential  equation;  GDL,  gas 
diffusion  layer;  PEM,  proton  exchange  membrane. 
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Fuel  cell  modeling  may  be  divided  into  two  categories— 
transport  modeling  and  system  modeling.  Transport  models  are 
used  when  specific  details  of  transport  phenomena  need  to  be 
studied.  They  are  comprehensive;  they  study  the  relevant  transport 
phenomena  occurring  within  the  fuel  cell,  and  solve  the  governing 
equations  in  the  multi-dimension  space-time  fuel  cell  domain.  The 
level  of  detail  and  complexity  often  determines  the  required  com¬ 
putational  resources.  The  convergence  time  of  a  transport  model 
may  range  from  minutes  to  days  depending  on  the  level  of  complex¬ 
ity  and  the  available  computational  resources,  which  may  range 
from  a  single  desktop/laptop  to  a  parallel  processing  network. 
Transport  modeling  has  been  used  to  conduct  parametric  studies 
of  fuel  cell  performance  [3],  and  to  investigate  design  issues  such 
as  fuel  composition,  water  and  thermal  management  [8],  catalyst 
micro-structure  [7],  the  effect  of  flow  field  geometries  [5,6],  and 
the  effect  of  using  different  types  of  membranes  [9-16]. 

System  models,  on  the  other  hand,  are  used  when  the  entire 
power  plant  containing  the  fuel  cell  stack  needs  to  be  analyzed. 
These  models  are  not  devoted  to  the  detailed  analysis  of  transport 
phenomena,  but  rather  the  fuel  cell  system  as  a  whole,  and  specifi¬ 
cally  how  the  fuel  cell  interacts  (sometimes  in  real  time)  with  other 
physical  and  virtual  components.  As  a  result  they  require  less  com¬ 
putational  resources.  System  models  have  been  used  in  real  time 
simulation  and  for  design  of  control  systems  [17],  and  thus  require 
fast  computations  being  performed  numerous  times.  Real  time 
models  require  these  computations  to  be  performed  up  to  1000 
times  every  second  [18].  Control  models,  employing  feedback  loops, 
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Nomenclature 

a  effective  surface  area  (m_1 ),  reaction  term  constant 

c  specific  heat  capacity  (J  kg-1  K_1 ),  constant  of  pro¬ 

portionality 

Di  j  diffusivity  of  gas  pair  i-j  (m2  s_1 ) 

F  Faraday  constant,  96487  C  mol-1 

i  current  density  (A  m-2 ) 

i0  exchange  current  density  (A  m-2 ) 

j  reaction  rate  (A  m-3 ) 

k  thermal  conductivity  (Wm-1K-1),  mathematical 

constant 

m  mole  fraction 

M  molar  mass  (kg  mol-1 ) 

n  reaction  parameter 

P  pressure  (Pa) 

R  reaction  term,  universal  gas  constant, 

8.3143  J  mol-1  K-1 
S  source,  entropy 

T  temperature  (K) 

u  velocity  (ms-1) 

Greek  letters 

a  charge  transfer  coefficient 

P  gas  permeability  (m2) 

y  concentration  parameter 

e  porosity 

fi  dynamic  viscosity  (Pas) 

p  density  (kg  m-3) 

o  electrical  or  ionic  conductivity  (S  m_1 ) 

(p  electrical  or  ionic  potential  (V) 

Subscripts  and  superscripts 
e  electrolyte  phase 

eff  effective 

ij  species  i,  j 

s  solid  phase 

T  thermal 


also  require  numerous  computations  to  be  performed  in  a  short 
space  of  time.  So  for  system  models,  it  is  necessary  to  greatly  sim¬ 
plify  the  problem  in  order  to  reduce  computational  requirements. 
Consequently,  system  models  are  usually  OD  (lumped  models)  or 
ID  [19].  Essentially  system  models  compromise  detail  for  speed. 
For  certain  applications,  especially  those  involving  control  systems, 
there  is  the  need  to  either  increase  the  level  of  detail  of  system  mod¬ 
els,  or  reduce  the  computational  requirements  of  transport  models. 
Mathematical  techniques  designed  to  reduce  the  computational 
requirements  of  transport  models  are  the  focus  of  this  paper. 

A  purely  mathematical  model  will  require  very  little  computer 
memory,  since  it  will  analytically  solve  the  fuel  cell  governing 
equations  resulting  in  a  closed  form  solution.  Flowever,  this  is  not 
realistic  due  to  the  highly  non-linear  nature  of  the  fuel  cell  gov¬ 
erning  differential  equations,  and  thus  closed  form  mathematical 
solutions  do  not  exist.  A  purely  computational  model  on  the  other 
hand,  does  require  a  high  amount  of  computer  memory.  Such  mod¬ 
els  use  finite  element,  finite  volume  or  finite  difference  methods  to 
approximate  the  real  problem.  Meshing  is  primarily  responsible  for 
the  large  memory  requirements.  Most  of  the  models  published  in 
the  open  literature  are  computational,  involving  very  little  if  any 
mathematical  solutions.  As  such,  very  few  analytical  models  have 
been  reported,  and  the  most  that  do  are  designed  to  be  stand-alone 
models,  not  incorporated  into  computational  models  [20-25].  Lit- 


Fig.  1.  PEM  fuel  cell  regions  and  computational  sub-domains. 


ster  and  Djilali  [26]  published  1 D  analytical  solutions  of  the  cathode 
catalyst  layer,  which  they  implemented  in  a  2D  model  of  an  air- 
breathing  fuel  cell.  Since  the  meshing  requirements  are  most 
demanding  in  the  catalyst  region  (see  Figs.  1  and  2),  an  analytical 
solution  of  the  catalyst  problem  incorporated  into  a  computational 
model  will  greatly  reduce  the  computational  requirements.  The 
present  work  considers  both  the  anode  and  cathode  catalyst  layers, 
and  implements  the  analytical  solutions  in  a  3D  fuel  cell  model. 

2.  Rationale 

Meshing  is  difficult  in  a  fuel  cell  for  two  reasons.  First,  there  are 
very  large  geometric  aspect  ratios,  since  the  dimensions  normal  to 
the  membrane  cross-section  (x  direction)  are  up  to  three  orders  of 
magnitude  smaller  than  the  other  dimensions.  Secondly,  even  in  the 
x  direction,  there  are  significant  scaling  discrepancies,  since  the  cat¬ 
alyst  layer  is  up  to  two  orders  of  magnitude  thinner  than  the  other 
regions.  So  meshing  becomes  very  difficult,  often  requiring  scaling 
to  create  the  optimum  balance  between  reducing  memory  require¬ 
ments  and  compromising  model  accuracy.  Fig.  1  shows  a  schematic 
of  the  fuel  cell  sub-domains.  This  figure  depicts  straight  channel 
flow  fields  showing  half  of  the  channel  cross-section  and  half  of 
the  ribs,  with  symmetry  being  assumed  across  they  extremes.  Fig.  2 
shows  the  same  sub-domains  as  Fig.  1,  but  with  the  finite  element 
mesh  displayed.  It  shows  how  dense  a  mesh  must  be  used  in  the 
catalyst  layers. 

The  catalyst  layer  is  a  very  thin  (~10  fxm)  three-dimensional 
region  sputtered  with  catalyst  particles,  and  with  electrochemical 
reaction  rates  varying  along  the  thickness  of  the  layer.  Modeling 
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the  CL  as  a  finite-sized  region  is  very  demanding  and  much  care 
must  be  taken  with  regards  to  computer  memory  requirements  and 
convergence.  Also,  depending  on  the  reaction  rate,  there  may  be  a 
very  high  rate  of  species  depletion  on  the  outer  surface  of  the  cata¬ 
lyst  region,  thus  confining  the  electrochemical  reactions  to  an  even 
thinner  layer  of  the  catalyst  region.  In  such  cases,  an  even  denser 
mesh  must  be  used  in  that  ultra  thin  region. 

Some  modelers  have  precluded  this  problem  by  treating  the 
catalyst  region  as  an  infinitesimally  small  interface,  rather  than 
a  finite-sized  region,  all  of  the  reaction  phenomena  being  reck¬ 
oned  on  this  interface.  Typically,  the  species  concentration  at  this 
interface  is  used  in  the  Butler-Volmer  equation  (or  Tafel  equation), 
for  which  the  effective  catalyst  surface  area  must  be  adjusted  to 
account  for  the  transformation  from  the  actual  3D  porous  catalyzed 
region  to  a  simplified  2D  catalyzed  interface.  One  problem  with  this 
approach  is  that  it  does  not  account  for  the  concentration  variation 
that  actually  exists  across  the  3D  catalyst  region. 

In  this  work,  we  develop  mathematical  solutions  to  the  phe¬ 
nomena  occurring  in  the  3D  catalyst  region.  These  solutions  are 
used  to  quantify  the  amount  of  species  consumption  or  gener¬ 
ation,  and  consequently  to  more  accurately  specify  the  flux  of 
species  occurring  at  the  GDL/CL  interface.  In  other  words,  these 
techniques  mathematically  transform  the  volumetric  catalyst  layer 
source  terms  into  interfacial  boundary  conditions  (flux  terms).  This 
treatment  precludes  the  problem  of  meshing  the  catalyst  layer,  thus 
reducing  the  computational  requirements  of  the  model,  while  at 
the  same  time  utilizing  a  more  accurate  3D  to  2D  transformation 
than  is  currently  practiced. 

3.  Model  equations 

The  fuel  cell  model  consists  of  a  highly  coupled  non-linear 
multi-physics  problem  involving  flow  in  porous  media,  diffu¬ 
sion  and  convection  of  species,  heat  transfer,  charge  transfer  and 
electrochemical  reactions.  These  phenomena  are  modeled  using 
the  Navier-Stokes  equations,  Darcy’s  law,  the  Stefan-Maxwell 
equations,  the  heat  equation,  Ohm’s  law,  and  the  Butler-Volmer 
equation.  These  governing  differential  equations  (GDEs)  are  listed 
below. 


r  Js. 

H2  “  ~2F 


where  j  represents  the  rate  of  proton  generation,  which  is  positive  at 
the  anode  and  negative  at  the  cathode.  Note  in  Eq.  (2),  the  term  P/RT 
is  used  for  the  total  molar  gas  concentration  to  avoid  an  obfuscation 
of  terminology,  and  is  assumed  to  be  locally  constant,  consistent 
with  the  earlier  simplifications,  m  is  the  mole  fraction,  and  S;  is  the 
molar  source  term.  The  following  source  terms  exist  only  in  the 
catalyst  layers,  where  electrochemical  reactions  take  place.  They 
quantify  the  consumption  of  oxygen  and  hydrogen,  as  well  as  the 
generation  of  water  and  heat,  all  as  a  function  of  the  reaction  rate, 
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(1) 

v  (i 

(2) 

(2) 

V  •  ( pcuT  -  pck  VT)  =  ST 

(3) 

V  •  is  =  V  •  (-CTseff  V0s)  =  -j 

(4) 

V  •  ie  =  V  •  (-<7eeff  V0e)  =  +j 

(5) 

The  boundary  conditions  pertinent  to  this  paper  are  the  conditions 
at  the  GDL/CL  and  the  CL/PEM  interfaces.  Species  enter  the  catalyst 
region,  react,  and  no  flow  occurs  across  the  membrane.  Across  the 
CL/PEM  boundary,  there  is  no  fluid  flow  and  no  diffusion  of  species. 
This  is  true  of  PBI  membranes.  The  Nation®  membrane  does  allow 
species  crossover,  especially  water.  Our  research  is  primarily  inter¬ 
ested  in  alternative  high  temperature  membranes,  such  as  PBI,  so 
our  results  focus  on  PBI  fuel  cells. 

4.  Analytical  techniques 

It  is  close  to  impossible  to  solve  the  full  fuel  cell  problem  using 
mathematical  techniques  alone  because  of  the  highly  non-linear 
nature  of  the  problem.  It  is  possible,  however,  to  solve  it  using 
numerical  techniques,  which  requires  extensive  meshing  and  com¬ 
putational  requirements,  depending  on  the  level  of  detail  desired. 
This  section  discusses  the  simplifications  of  the  physical  prob¬ 
lem  necessary  for  a  purely  mathematical  treatment  of  the  catalyst 
layer  problem,  and  the  development  of  those  approaches.  These 
mathematical  solutions  are  intended  to  replace  the  need  for  com¬ 
putational  modeling  of  the  catalyst  layer.  The  volumetric  CL  source 
terms  are  transformed  into  interfacial  boundary  conditions  to  be 
specified  at  the  GDL/CL  interfaces.  The  end  result  will  be  a  model 
that  employs  mathematical  rather  than  computational  solutions 
of  the  catalyst  layer  problem.  This  will  be  called  a  semi-analytical 
model  and  is  shown  schematically  in  Fig.  3. 

When  this  model  is  implemented  in  a  computational  model,  the 
catalyst  layer  sub-domains  are  merged  with  the  membrane  to  form 
one  larger  sub-domain  where  all  ionic  transport  takes  place.  This 
will  minimize  the  error  in  determining  the  effective  ionic  resistance 
of  the  cell. 

The  following  simplifications  are  made  in  order  to  enable  math¬ 
ematical  solutions  to  the  catalyst  layer  equations. 


diffusion,  we  neglect  the  convection  terms  in  the  catalyst  layer 
GDEs. 

We  assume  anisotropic  behavior  in  the  catalyst  region.  Since  the 
catalyst  region  is  thin,  we  assume  that  the  transport  of  species 
occurs  primarily  in  the  x  direction  (normal  to  the  cross-section). 
Therefore,  we  neglect  transport  in  they  and  z  directions,  which 
enables  a  pseudo  ID  treatment  of  the  catalyst  problem.  Note 
that  this  restriction  applies  only  in  the  catalyst  layers  and  not 
the  gas  diffusion  layers  or  the  gas  channels. 

(3)  We  further  assume  that  since  the  catalyst  region  is  thin  in 
the  x  direction,  changes  in  temperature,  pressure  and  activa¬ 
tion  overpotential  occurring  across  the  layer  can  be  neglected. 
Essentially,  we  assume  these  quantities  only  vary  in  the  y  and 
z  directions  in  the  catalyst  layer,  and  there  is  no  interaction 
between  state  properties  at  different  (y,  z)  locations,  consis¬ 
tent  with  assumption  (2).  Again,  these  restrictions  apply  only 
to  the  catalyst  layers.  It  should  be  noted  that  we  do  not  assume 
isothermal  conditions,  only  constant  x  direction  temperature 
along  the  thin  catalyst  layer  at  each  (y,  z)  location. 

(4)  Species  concentration  and  reaction  rates  are  the  only  param¬ 
eters  that  vary  significantly  across  the  catalyst  layer  in  the  x 
direction.  We  do  not  assume  constant  concentration  of  species 
or  constant  electrochemical  rates  across  the  catalyst  layer,  since 
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Fig.  3.  Schematic  and  scope  of  the  mathematical  techniques  presented. 


assumptions  of  this  nature  are  far  from  realistic.  Instead  we  offer 
analytical  solutions  of  the  x  variation  in  concentration  at  each 
(y,  z)  location. 

With  these  assumptions,  the  Butler-Volmer  Eq.  (6)  reduces  to 

j(x,y,z)  =  a(y,z)m(x,y,z)n  (11) 

where  the  species  subscript,  i,  is  dropped,  and  the  value  a  is 
an  x  direction  constant  representing  all  the  other  terms  in  the 
Bulter-Volmer  equation.  From  this  point,  the  (x,  y,  z)  notation  will 
be  dropped.  The  typical  values  for  n  in  the  PEM  fuel  cell  are  1 
(cathode)  and  1/2  (anode). 

We  see  that  the  volumetric  reaction  rate  depends  solely  on 
the  species  concentration— hydrogen  at  the  anode,  oxygen  at  the 
cathode.  Further,  the  species  consumption  source  terms  are  linear 
functions  of  j. 

5  =  -kj  =  -kamn  (12) 

where  k  is  a  positive  constant. 

The  species  conservation  Eq.  (2)  in  the  catalyst  region  becomes 
AuVm  =  0=^DV2m  +  S 

K1  K1 


(13) 


d^m  /fcom  2 

dx2  V  PD  J 


offer  approximate  analytical  solutions  based  on  typical  conditions 
that  often  occur  in  PEM  fuel  cells.  The  first  case  is  an  approxima¬ 
tion  for  fast  reaction  rates.  The  second  case  involves  a  Taylor  series 
approximation  to  translate  the  power  term  to  a  linear  term.  The 
third  case  applies  to  very  slow  reaction  rates,  where  the  species  con¬ 
centration  does  not  change  significantly  across  the  catalyst  layer. 
These  cases  will  then  be  compared  for  accuracy. 

4  A.  Case  l —fast  reaction  rates 

When  the  species  consumption  rate  is  very  high,  the  concentra¬ 
tion  quickly  decreases  to  zero  within  the  catalyst  layer.  This  may 
happen  at  high  current  densities,  or  when  a  very  active  catalyst 
layer  is  present.  Paradoxically,  in  such  cases,  the  thin  catalyst  layer 
can  be  treated  as  infinitely  large  compared  to  the  length  scale  of 
the  electrochemical  reaction. 

Using  the  substitution, 


nr 


(14) 


dx 

Eq.  (14)  becomes, 

— —  b2  mn 
if3  dm 

This  can  be  solved  using  separation  of  variables  to  obtain 


(17) 


(18) 


Eq.  (14)  is  the  simplified  x  direction  species  transport  equation, 
which  is  solved  at  each  (y,  z)  location.  This  equation  is  defined  in  the 
catalyst  layer  (0  <x<  r),  where  r  is  the  catalyst  layer  thickness,  and 
x  is  measured  from  the  GDL  to  the  CL.  Note  that  for  the  cathode  CL, 
this  co-ordinate  is  directed  in  the  opposite  direction  to  that  shown 
in  Figs.  1  and  2. 

The  concentration  at  x  =  0  is  known  from  the  computational 
solution  of  the  GDL  problem,  while  insulation  conditions  exist  at 
x  =  r.  The  mathematical  representations  of  these  boundary  condi¬ 
tions  are 

m(0,y,z)  =  m0(y,z)  (15) 

A(m(z,y,z)  =  0  (16) 

For  the  case  of  n  =  1,  the  solution  to  Eq.  (14)  is  very  straightforward. 
For  n  =  1/2,  analytical  closed  form  solutions  do  not  exist.  Instead,  we 


1  dm  ^(2  m1+nV^2 

i[r  ~  dx  ~  ll+ny 


(19) 


The  constant  of  integration  is  set  to  zero  in  order  to  satisfy  the 
condition  that  both  the  concentration  and  flux  are  zero  when  x  =  r. 
Since  from  Fields  law,  diffusion  is  proportional  to  the  negative  of  the 
first  derivative  of  concentration,  we  use  the  negative  square  root  in 
Eq.  (19).  From  Eq.  (19),  we  can  obtain  the  boundary  flux  in  terms  of 
the  mole  fraction  at  x  =  0. 


dm 


dx 


x=0 


-b 


2  mj+n 

TTn“ 


1/2 


(20) 


PD  dm  (  2kaPDrrig+n  \  1/2 

RT  dx  x=0\  RT(l  +  n)  J 


(21) 
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This  term  represents  the  flux  of  species  leaving  the  GDL  and  enter¬ 
ing  the  CL  at  position  (y,  z)  in  terms  of  the  species  mole  fraction  at 
(y,  z).  This  term  now  replaces  the  volumetric  generation  term  in  the 
catalyst  region.  Instead  of  an  insulation  condition  being  specified 
at  the  CL/PEM  interface,  a  flux  (or  convection)  condition  is  speci¬ 
fied  at  the  GDL/CL  interface,  which  accounts  for  the  electrochemical 
reactions  occurring  in  the  CL. 

The  solution  can  be  further  developed  to  get 


mW2_ml-n/2 


a  1  n 

2V2(1  +n)1/2* 


(22) 


This  is  valid  for  0  <x<  r.  The  mole  fraction,  m  falls  to  zero  when 


[2(1  +n)]1/2  n/2 

a(l  -  n)  0 


(23) 


We  expect  the  case  1  solution  to  provide  accurate  results  when  the 
catalyst  layer  is  larger  than  the  RHS  of  Eq.  (23),  which  provides  an 
estimate  of  the  thickness  of  the  catalyst  layer  actually  being  utilized 
for  electrochemical  reactions. 

All  of  the  other  source  terms  can  be  written  in  terms  of  the 
species  source  term.  Suppose  a  particular  source  term  is  linearly 
related  to  the  species  source  term,  such  that 


S  =  cb2mn 


(24) 


where  c  is  the  constant  of  proportionality.  To  find  the  total  genera¬ 
tion  (or  consumption)  of  that  particular  source  term,  S  needs  to  be 
integrated  over  the  catalyst  region.  The  result  would  be, 


cb2  f 

m 

n  dx  , 
n  —  dm 
dm 

J  m= 

--m0 

*r 

0)n 

/  1  +  n  \1/: 
\2in1+n/ 

J  CD=COq 

dm 

C  dx 

x=0 

cb(  . 

2  \ 

1/2 

mn+n)/2 

VI  +nj 

dm 
— c  — 

dx 

x=0 

dm 


(25) 

(26) 

(27) 

(28) 
(29) 


Not  surprisingly,  this  is  the  same  result  that  was  previously 
obtained  in  Eq.  (21).  So  all  source  terms  (which  are  proportional 
to  the  species  source  term)  can  be  transformed  and  written  as  lin¬ 
ear  functions  of  the  species  boundary  flux  with  the  same  constants 
of  proportionality. 


4.2.  Case  2— Taylor  series  approximations 


In  cases  where  the  reactions  are  not  fast  enough  to  totally  con¬ 
sume  the  reactants,  the  case  1  solution  would  lose  accuracy.  In  fact, 
we  would  expect  it  to  overestimate  the  reaction  rate.  In  such  cases, 
we  may  opt  for  an  “exact”  solution.  If  n  =  1,  the  differential  equation 
is  easily  solvable.  For  0<n<  1,  we  can  use  the  binomial  expansion 
to  linearize  the  power  term,  resulting  in  a  solvable  equation. 

mn  =  [  1  -  (1  -  m)]n  ^  1  -  n  (1  -  m)  =  nm  +  (1  -  n)  (30) 


This  approximation  obviously  works  best  when  the  mole  fraction 
is  close  to  unity,  and  is  actually  exact  when  n  =  1.  The  solution  is 

dm  _  ^^2  1  -  exp(2brn1/2)  /nm0  +  1  -  n\ 
dx  x=0  n  1  +  exp(2brn1/2)  V  n  ) 


b  1 -exp(2brn1/2) 
n1/2  1  +  exp(2brn1/2)  0 


(31) 


PD  dm  _  /kaPD\ 1  -  exp(2farn1''2)  „ 

RTd7x=0  =\1W)  1  +  exp(2farn1/2)m° 


(32) 


4.3.  Case  3— slow  reaction  rates 


If  the  reaction  rates  are  very  slow,  the  concentration  of  the 
species  does  not  vary  significantly  across  the  catalyst  layer.  This 
is  sometimes  the  case  with  hydrogen.  Since  the  reaction  rate  con¬ 
stant  for  the  oxidation  of  hydrogen  is  significantly  faster  than  the 
rate  of  oxygen  reduction,  the  latter  is  the  rate  controlling  step.  So 
although  the  hydrogen  is  readily  available  for  reaction,  hydrogen 
only  reacts  at  a  rate  necessary  to  provide  a  sufficient  supply  of  pro¬ 
tons  for  the  oxygen  reaction.  This  rate  is  much  less  than  the  rate  at 
which  hydrogen  is  capable  of  reacting,  and  as  a  result  the  hydrogen 
concentration  does  not  vary  significantly  in  the  x  direction  across 
the  catalyst  layer.  In  such  cases,  we  may  assume  that  the  concen¬ 
tration  is  constant  across  the  catalyst  layer.  The  amount  of  species 
that  reacts,  which  is  the  same  as  the  amount  that  diffuses  across 
the  interface,  is  simply 


PD  dm 
RT  dx 


=  -karm^ 

x=0 


(33) 


4.4.  Summary  of  all  cases 

At  the  anode,  n  =  1/2,  and  at  the  cathode,  n  =  1. 


(aAE  \ 

(  acF  \ 

exp  (rH 

-exp  (-wnj 

j^anode  _ 


1 

2 F 


^cathode  _ 

-  4  F 


(34) 

(35) 

(36) 


The  following  quantities  represent  the  boundary  condition  at  the 
GDL/CL  interface  that  replace  the  volumetric  CL  source  terms. 


Case  1 : 

(  4kaPDH,  \  1/2  3/4 

flUX(H2)=  < 

flux(0  2)=(^y/2m02 


(37) 


(38) 


Case  2: 


flux(H2)  = 


1  -exp[r(2fcaRT/PDH2)1/2]  (2kaPD^2m^ 
1  +  exp[r(2 kaRT/PDH2  )1/2 ]  \  RT 


1/2 


(39) 


flux(02)  = 


1  -  exp[2r (kaRT/PDo2)y2]  f  kaPD02 
1  +  exp[2T(kaRT/PDo2 )1/2]  V  RT 


1/2 

m02 

(40) 
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Fig.  4.  Concentration  profile  for  ID  validation  (a  =  10-3,  n  =  l). 


Case  3: 

flux(H2)  =  karm^ 

(41) 

flux(02)  =  kaxmo2 

(42) 

The  remaining  boundary  conditions  can  be  written  as  linear  factors 
of  the  hydrogen  and  oxygen  fluxes.  ia  and  ic  are  respectively,  the 
electronic-to-ionic  current  transfer  at  the  anode  and  the  ionic-to- 
electronic  current  transfer  at  the  cathode,  q  is  the  heat  of  reaction 
evaluated  at  the  cathode  GDL/CL  interface. 

flux(H20)  =  -2flux(02) 

(43) 

ia  =  2Fflux(H2) 

(44) 

ic  =  4Fflux(02) 

(45) 

(  TAS\  . 

Qrxn  =  yjl  —  np  J  <c 

(46) 

5.  Results  and  discussion 


will  be  applied  to  a  full  3D  model  of  a  PEM  fuel  cell  equipped  with 
a  PBI  membrane. 


5.1.  Validation  with  a  simplified  problem 


Consider  a  simplified  conduction-convection  problem  gov¬ 
erned  by  species  conservation  and  Darcy’s  law.  All  units  are 
consistent  with  the  SI  system,  although  not  very  important  in  this 
simplified  example. 


-10 


-6  d2C 


dC  D 

,  9  +  M-J—  —  R 
dx2  dx 


^=_10-4^=R 

dx  dx2 


where, 


/O,  -1  <x<  0  1 

|-aCn,  0  <  x  <  0.2  J 

The  boundary  conditions  are 


(47) 

(48) 


(49) 


The  formulations  derived  in  Section  4  are  now  tested 
for  validation.  First  they  will  be  applied  to  a  simple  ID 
convection-conduction  problem  with  a  source  term.  Secondly,  they 


C(-l)=l 

(50) 

^5 

1 

II 

o 

U1 

(51) 

Fig.  5.  Concentration  profile  for  ID  validation  (a  =  10-5,  n  =  1). 
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Fig.  6.  Concentration  profile  for  ID  validation  (a  =  10-3,  n  =  1/2). 


AC(0.2)  =  0  (52) 

^P(0.2)  =  0  =  U(0.2)  (53) 

This  is  first  solved  using  a  finite  element  approach  with  the  reaction 
region  (0<x<0.2)  meshed  and  the  spatial  source  term,  R  applied. 
This  full  computational  solution  is  called  C.  It  is  then  solved  using 
the  analytical  catalyst  layer  treatment  described  in  Section  4.  These 
solutions,  presented  for  cases  1,  2  and  3,  are  respectively  called 
Cl,  C2  and  C3.  The  two  quantities  that  are  varied  are  a  and  n.  a 
represents  the  rate  of  reaction  term— a  high  value  representing  a 
fast  reaction  rate,  and  n  represents  the  order  of  the  reaction. 

Consider  the  case  when  n  =  l  and  a  =  10-3  (Fig.  4).  This  repre¬ 
sents  a  fast  reaction  rate— evident  by  the  fact  that  the  concentration 
C  falls  to  zero  within  the  catalyst  region  (0<x<0.2).  C2  matches 
almost  exactly  with  C,  while  Cl  and  C3  slightly  underestimate  the 
concentration.  With  n  =  1,  C2  provides  an  exact  solution  to  the  cata¬ 
lyst  layer  equations.  Therefore,  it  is  not  surprising  that  it  accurately 
simulates  C. 

Setting  a  =  10-5,  which  represents  a  slower  reaction  rate,  results 
in  an  exaggeration  of  this  result.  Again,  and  expectedly  as  Fig.  5 
shows,  C2  matches  almost  exactly  with  C.  Now  Cl  and  C3  underesti¬ 
mate  the  concentration  by  an  even  greater  amount.  This  is  expected 
since  Cl  works  best  when  the  reaction  rate  is  fast.  Thus,  the  higher 
the  value  of  a,  the  better  is  the  accuracy  of  Cl.  In  both  cases,  C3 
has  a  small  error.  C3  assumes  a  constant  concentration  across  the 
catalyst  layer  which  is  not  the  case  in  this  problem.  As  a  result, 
it  must  mathematically  reduce  the  concentration  everywhere  to 
compensate  for  this  inaccuracy.  Most  PEM  fuel  cell  models,  which 
use  an  interfacial  catalyst  layer  treatment,  employ  a  method  which 
more  closely  resembles  C3  than  any  of  the  other  two.  Note  that  this 
formulation  does  not  always  provide  accurate  results. 

Now  setting  n  =  1/2  and  a  =  10-3  (fast  reaction  rate),  Fig.  6  shows 
that  Cl  provides  the  best  estimate  of  C.  C2  is  only  accurate  when  C  is 
close  to  1,  which  is  not  the  case  here.  So  C2  slightly  underestimates 
the  concentration.  For  similar  reasons  as  above,  C3  also  underes¬ 
timates  the  concentration.  Therefore,  it  is  not  surprising  that  Cl 
provides  such  accurate  results,  since  it  assumes  very  fast  reaction 
rates. 

Setting  a  =  10-6  (Fig.  7),  which  represents  a  slow  reaction  rate, 
Cl  now  underestimates  the  concentration  since  it  imposes  a  higher 
reaction  rate  than  is  warranted.  In  this  case,  C2  and  C3  provide 
accurate  results.  C2  is  accurate  because  the  concentration  is  close 


to  1,  thus  rendering  the  Taylor  series  approximation  valid.  C3  is 
accurate  because  the  actual  concentration,  C  in  the  catalyst  region 
does  not  change  significantly.  For  the  anode  reaction  in  PEM  fuel 
cells,  this  is  the  most  realistic  scheme.  So  from  the  simplified  valida¬ 
tion  test,  it  seems  that  C2  and  C3  would  provide  the  most  accurate 
scheme  for  the  anode  reaction,  while  C2  appears  to  be  the  best 
choice  at  the  cathode  where  n  =  1. 

5.2.  Validation  with  full  3D  PEM  fuel  cell  model 

These  analytical  formulations  are  now  incorporated  in  our  pre¬ 
vious  3D  computational  PBI  fuel  cell  model  [12],  and  results  are 
compared.  The  previous  model  is  fully  computational,  while  the 
new  implementation  is  semi-analytical.  All  numerical  constants 
and  operating  conditions  remain  unaltered.  These  constants,  listed 
in  [12],  are  not  repeated  here  since  they  are  not  relevant  to  the 
findings  of  this  paper.  It  is  only  the  comparisons  of  the  various  for¬ 
mulations  with  the  full  computational  model  that  is  of  essence. 
All  concentration  profiles  are  shown  for  an  operating  cell  voltage 
of  0.3  V.  The  notation  HiOj  means  that  the  case  i  formulation  is 
used  on  the  hydrogen  side  (anode),  while  the  case  j  formulation 
is  used  at  the  oxygen  side  (cathode).  H0O0  refers  to  the  original 
computational  model. 

Fig.  8  shows  the  IV  curve  of  the  original  computational  model 
and  five  of  the  semi-analytical  solutions.  FI302,  FI303  and  FI202 
all  provide  a  very  close  match  with  FI0O0,  such  that  they  are 
indistinguishable  from  each  other  on  the  graph.  Both  HI 02  and 
H301  produce  inaccurate  results  since  they  require  that  the  con¬ 
centration  of  hydrogen  and  oxygen,  respectively,  fall  to  zero  in 
the  catalyst  layer.  As  a  result,  they  both  impose  higher  reaction 
rates  than  are  warranted  by  the  actual  problem.  H301  is  interest¬ 
ing.  At  low  current  densities  (or  low  reaction  rates),  the  deviation 
from  HOOO  is  very  large.  However,  as  the  current  density  increases, 
the  error  becomes  smaller.  This  is  expected  since  at  higher  cur¬ 
rent  densities,  the  depletion  of  oxygen  becomes  more  pronounced. 
It  can  be  expected  that  in  cases  where  the  oxygen  concentra¬ 
tion  is  fully  depleted,  the  01  formulation  would  provide  accurate 
results. 

Fig.  9  shows  the  hydrogen  mole  fraction  along  the  gas  channels 
and  gas  diffusion  layers  for  the  full  computational  model,  HOOO. 
The  mole  fraction  in  the  anode  catalyst  layer  is  not  shown  for  the 
purpose  of  comparisons.  The  mole  fraction  varies  from  an  inlet  of 
0.963  to  a  lowest  value  of  0.961.  In  this  case,  the  hydrogen  concen¬ 
tration  does  not  show  much  variation,  neither  in  the  z  direction 


(A)  A 


D.F.  Cheddie,  N.D.H.  Munroe  /  Journal  of  Power  Sources  183  (2008)  164-173 


171 


x 


Fig.  7.  Concentration  profile  for  ID  validation  (a  =  10-6,  n  =  1/2). 


Fig.  8.  Polarization  curves  (HOOO,  H102,  H202,  H302,  H301,  H303). 
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Fig.  9.  Hydrogen  mole  fraction  (HOOO). 
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Fig.  10.  Hydrogen  mole  fraction  (H202,  H203,  H303). 
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nor  in  the  x  direction.  H202,  H302  and  H303  all  produce  results 
which  are  barely  distinguishable  from  each  other.  Consequently, 
only  one  diagram  (Fig.  10)  is  used  to  display  all  of  them.  These 
solutions  predict  a  similar  concentration  profile,  but  with  a  low 
mole  fraction  of  0.956,  which  is  a  slight  deviation  from  the  0.961  of 
the  computational  model.  However,  in  the  context  of  the  hydrogen 
concentration  that  does  not  vary  significantly,  this  deviation  is  not 
considerable. 

Fig.  11  shows  the  oxygen  mole  fraction  for  H0O0,  with  the  mass 
fraction  ranging  from  0.126  to  0.194.  This  represents  a  much  more 
significant  concentration  variation,  as  is  expected  with  oxygen. 
Fig.  12  shows  the  equivalent  plot  for  H202,  H302  and  H303,  respec¬ 
tively,  which  as  was  the  case  with  hydrogen,  are  indistinguishable 
from  each  other.  The  analytical  solutions  match  the  oxygen  mole 
fraction  almost  perfectly,  with  the  computational  solution,  showing 
virtually  no  difference  in  range  or  variation. 

These  results  show  that  the  mathematical  formulations  all  pro¬ 
vide  a  very  accurate  representation  of  the  full  computational  model, 
but  it  greatly  reduces  the  computational  memory  requirements 


and  computation  speed,  since  there  is  no  need  for  meshing  of  the 
catalyst  layer.  The  full  computational  model  requires  5350  second- 
order  tetrahedral  elements,  while  the  semi-analytical  models  only 
require  1910  finite  elements.  The  computational  model  requires 
24  min  to  gradually  converge  upon  a  solution,  while  the  semi- 
analytical  models  require  6  min,  all  starting  with  identical  initial 
estimates.  All  computations  were  performed  on  a  Windows  XP  plat¬ 
form,  with  a  2  GHz  Intel  Pentium  4  processor  and  512  MB  RAM, 
using  the  multi-physics  software  FEMLAB  3.1  i®. 

This  becomes  particularly  useful  when  dealing  with  larger  scale 
models.  The  present  model  only  simulated  one-half  channel  and 
half  rib,  and  assumed  y  direction  symmetry.  If  it  were  required 
to  simulate  serpentine  flow  fields,  for  instance,  the  computational 
requirements  would  be  much  greater.  The  formulations  presented 
here  can  help  to  reduce  that  requirement  and  allow  larger  problems 
to  be  modeled  with  less  computational  resources.  The  techniques 
presented  here  can  also  be  applied  to  2D  transport  models,  and 
may  simplify  them  enough  to  be  used  in  system  models  involving 
control  systems  and  other  applications. 


Fig.  11.  Oxygen  mole  fraction  (H0O0). 
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Fig.  12.  Oxygen  mole  fraction  (H202,  H203,  H303). 


6.  Conclusions 

Analytical  solutions  were  presented  for  the  transport  and  elec¬ 
trochemical  problem  in  PEM  fuel  cell  catalyst  layers  at  both  the 
anode  and  cathode.  These  solutions  transformed  the  volumetric 
reaction  terms  to  boundary  fluxes,  thereby  eliminating  the  need  to 
mesh  the  catalyst  layer.  The  result  was  a  model  that  required  less 
computational  memory  and  converged  in  less  time. 

These  semi-analytical  solutions  matched  very  well  with  a  full 
computational  model  in  terms  of  the  polarization  results,  hydrogen 
and  oxygen  concentration.  The  results  showed  that  the  analytical 
techniques  did  not  compromise  the  accuracy  of  the  model.  It  was 
further  suggested  that  these  techniques  could  reduce  the  compu¬ 
tational  requirements  of  2D  transport  models,  whereby  allowing 
them  to  be  incorporated  into  system  models. 

The  techniques  presented  in  this  paper  applied  primarily  to 
high  temperature  PEM  fuel  cells,  using  alternative  membranes  to 
Nation®,  e.g.  PBI. 
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